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The aging part R ag (t, s) of the impulsive response function of the two dimensional ferromagnetic 
Ising model, quenched below the critical point, is studied numerically employing a new algorithm 
without the imposition of the external field. We find that the simple scaling form R ag (t,s) = 
s _ ( 1 + a ) f(t/s), which is usually believed to hold in the aging regime, is not obeyed. We analyse the 
data assuming the existence of a correction to scaling. We find a = 0.273 ± 0.006, in agreement 
with previous numerical results obtained from the zero field cooled magnetization. We investigate in 
detail also the scaling function f(t/s) and we compare the results with the predictions of analytical 
theories. We make an ansatz for the correction to scaling, deriving an analytical expression for 
R ag (t,s). This gives a satisfactory qualitative agreement with the numerical data for R ag (t,s) and 
for the integrated response functions. With the analytical model we explore the overall behavior, 
extrapolating beyond the time regime accessible with the simulations. We explain why the data for 
the zero field cooled susceptibility are not too sensitive to the existence of the correction to scaling 
in R ag (t,s), making this quantity the most convenient for the study of the asymptotic scaling 
properties. 

Jcorberi@sa.infn.it flippicllo@sa.infn.it §zannetti@na. infn.it 
PACS: 05.70.Ln, 75.40.Gb, 05.40.-a 

I. INTRODUCTION 

A powerful tool in the study of slowly relaxing systems is the extension of the fluctuation dissipation theorem (FDT) 
to off-equilibrium conditions. The relation between the response and the autocorrelation function has been shown 
to encode basic properties of the dynamics Q and of the equilibrium state 0- However, while the autocorrelation 
function C(t,s) has been studied since a long time and, in some systems, is well understood, the knowledge of the 
response function R(t, s) remains comparatively poor. 

In the prototypical case, considered here, of non disordered coarsening systems quenched below the critical tem- 
perature, general properties of C(t, s) and R(t, s) may be inferred from the structure of the configurations. It is well 
known, in fact, that in the late stage of phase ordering the interior of the growing domains is equilibrated, while their 
boundary is out of equilibrium. Accordingly, a distinction between bulk and interface fluctuations can be made. For 
the autocorrelation function this leads to the splitting 0. [J 0- E| 

C(t,s) = C st (t-s) + C ag {t,s) (1) 

where C s t(t— s) is the contribution from the interior of domains, behaving as the equilibrium autocorrelation function 
C eq (t — s) of the ordered state at the final temperature T of the quench, while C ag (t, s) is the off-equilibrium aging 
contribution, coming from the interfaces. Similarly, for the response function one has 

R(t,s) = R st (t-s) + R ag (t,s) (2) 

and the stationary contributions R s t(t — s) and C s t(t — s) are related by the FDT 

TR at (t-s)= dC ^- a) . (3) 

For s large, C ag {t, s) obeys the simple aging scaling form 

C ag {t,s) = s- b g{x) (4) 

with b = and x = t/s [HEUl] Taking s = and t large, C ag (t, 0) decays algebraically as C ag (t, 0) ~ t~ x l z , where 
A is the Fisher-Huse exponent and z is the dynamic exponent regulating t he g rowth of the characteristic size 
Lit) ~ t 1 / 2 of domains. Assuming that the same behavior holds also for s > [l(| gives g(x) ~ x~ x / z , for x ^> 1. 
These properties of C(t,s), as mentioned above, are well understood and well documented both by analytical and 
numerical results. 
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By analogy, the response function is also expected to obey, for large s, a simple aging scaling form 



Rag(t,s) = 8-^+4 f(x) 



(5) 



with 




investigating R(t, s) is considerably more difficult than C(t, s). The reason is that R(t,s) is much more noisy than 
C(t,s). A possible way around this difficulty is to study the less noisy integrated response functions (IRF), such as 
the zero field cooled magnetization (ZFC) 



However, if resorting to the IRF does indeed cut down the noise, it has turned out that recovering R(t, s) from the IRF 
is a delicate task which involves more than one subtlety (see Ref. an d Section[n]of this paper). As a result, so far, 
in the literature no consensus has been reached neither on the exponent a nor on the form of the scaling function f(x). 
Therefore, it seems that the issue of the properties of R(t, s) can be settled o nly by goin g to its direct measurement. 
Recently, this has become feasible after the introduction of new algorithms [TsL GjjTSll f° r the computation of the 
response function without applying the external perturbation. This speeds the simulation to such an extent that the 
direct measurement of R(t, s) has become accessible [Til lljL I20I l2lj] . 

In this paper we present results for R(t, s) in the two dimensional Ising model obtained with our algorithm |20j| . 
The outcome is quite rich and interesting. We find the unexpected result that in the d = 2 Ising model the simple 
scaling form (J5J is not enough to represent the data, but a large correction term is needed, up to the longest times we 
have reached in the simulations. The quality of the data allows to detect and analyze this correction and, by taking 
it properly into account, to make statements on the exponent a and the scaling function f(x) appearing in Eq. 

The paper is organized as follows: In Section [H] we summarise the problem of the exponent a. In Section ITTT1 we 
recall the basic features of the measurement of the response function without applying the external perturbation and 
the method we use to isolate R ag (t,s) in Eq. ©. In Section |W] we present and discuss the numerical results for 
R ag (t, s). In particular, we show the existence of a large preasymptotic scaling correction, we extract the value of the 
exponent a and we analyze the scaling function f(x). In Section we investigate the implications of the structure 
of R ag (t,s) on the properties of the IRF and we compare our predictions with the numerical simulations of these 
quantities. This gives further insight into the problem of the retrieval of R(t,s) from the IRF. Finally, in Section IVT1 
we draw the conclusions. 



Before entering the presentation of data and results, we recall briefly what is the problem with the exponent a in 
Eq. ©. 

• The straightforward substitution of Eq. JSJ into Eqs. Q and JSJ) yields the scaling forms of the aging parts of 




(7) 



or the thermoremanent magnetization (TRM) 




(8) 



II. THE EXPONENT a 



the IRF 



Xag(t,t w ) — t w x F x (y) 



(9) 



and 



Pag(t,t w ) =t w a "F p {y) 



(10) 



with 



a x = a p = a. 



(11) 



and y = t/t 
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• an intuitively appealing argument, originally introduced by Barrat [23, predicts that Xag(t,t w ) ought to be 
proportional to the density of defects at the time t w , which goes as L(t w )~ n , with n = 1 and n — 2 for scalar 
and vectorial order parameter, respectively. Using Eq. pill , this leads to the dimensionality independent result 

a x = a p = a = n/z (12) 

since, we recall, the exponent z in the quenches below Tc does not depend on d. For systems without conservation 
of the order parameter, as considered in this paper, z = 2 p|. 

• contrary to the previous statement, all available analytical results, which include the exact solutions of the large 
N model £| and of the d = 1 Ising model 0, 0] , as well as approximate calculations flU ITU IT(| , show that a 
depends linearly on dimensionality according to 

n d-d L 
z d ( j - d L 

where rf^ is the lower critical dimension of the statics and djj > 4 is a dimensionality whose significance 
becomes clear as soon as ZFC is considered. In fact, from these analytical calculations it comes out that the 
relation between a x and a is not as simple as in Eq. 1)11(1 . Rather, it must be replaced by 

as in Eq. |^3|) for d < du 

with logarithmic corrections for d — djj (14) 
for d> djj. 

This behavior of a x is due to the existence of an irrelevant variable in Xag(t,t w ), which becomes dangerous 
for d > du ln the large N model the above formula holds with (n = 2,d^ — 2,du = 4), while in the 

approximate calculations with scalar order parameter of Refs. 0,0,^1 ( n — E = ^-,du = 2). 

• The problem of the exponent a, then, is whether Eq. ((13(1 is a peculiarity of just those cases where analytical 
results are available, the generic behavior being that of Eq. I(12|) . or, viceversa, it is Eq. 1(13(1 that captures the 
generic behavior, revealing thereafter that the argument leading to Eq. ((12|) does to miss some important feature 
in the mechanism of the response. 

In order to answer the question one has to investigate as many systems as possible, with the aim of putting together 
the generic picture. This can be done only by numerical methods. We have carried out such a program performing 
simulations and measuring a x in several systems with conserved and non conserved dynamics, both with scalar and 
vectorial order parameter, at different dimensionalities j^. The large body of results that we have obtained does 
indeed indicate, quite convincingly in our opinion, that Eq. ((14() is of general validity, with djj = 3 |24| and du = 4 
for systems with scalar and vector order parameter, respectively. 

This conclusion is challanged, mainly on the basis of the measurement of a p from TRM in the Ising model with 
d = 2 and d = 3 [2{|, which seems to agree with Eq. 1(12(1 . In other words, the investigations of ZFC and TRM seem 
to produce different results. We have explained in detail elsewhere [T^ | that the difference is only apparent and is 
due to the fact that while the data for ZFC are asymptotic, those for TRM are not. To us, the interesting question 
remaining open is not anymore what is the value of a, but what is the physics behind Eq. ((13(1 . However, as of yet 
this is not a shared conclusion prl [271. l28l ] . So, what emerges from this brief account of the problem is that working 
with the IRF is a hairy business, because it has opened the somewhat intricate problem of why ZFC and TRM seem 
to give different results. Therefore, in order to make progress, a fresh starting is needed. In this respect, the direct 
numerical study of R(t, s), with the new algorithms, seems well suited to the task. We shall concentrate on the d = 2 
Ising model, where discrimination between Eqs. 1(12(1 and ((1311 ought to be easier, since the difference in the predicted 
values of a is quite large 

1/2 from Eq. JT2J) ^ 
1/4 from Eq. ^ ' 

III. THE ALGORITHM 
A. Measurement of the response without applying a perturbation 

The basic idea, in the new methods [l8L Il9l [20] for the measurement of the response function without applying 
the external perturbation, is to relate R(t, s) to some correlation function of the unperturbed dynamics, much in the 
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same way as in the equilibrium FDT. In this paper we will use our method introduced in |20| . since we have checked 
that it is numerically more efficient. Let us briefly describe it, referring to [2(| for details. 

We consider a spin system with Hamiltonian W[er], where [a] is a generic configuration of the spin variables Oi = ±1. 
The response function is defined by 



. . _ 1 dladt)) 

R(t, s = lim i-^i 

v ' a s ^o As dhi 



(16) 



h=0 



where hi is a magnetic field acting on the i-th site during the time interval [s, s + As], and the right hand side does not 
depend on i due to space translational invariance. Computing (crj(t)) by means of the master equation and inserting 
the result into Eq. (|16|l . one arrives at 



TR(t, s) = - lim 

y 1 J 2 As^O 



C(t,s + As)-C(t,s) 
As 



{ai(t- As)Bi(s)) 



where Bi enters the evolution of the magnetization according to 

d(ai(t)) 



dt 



(Bi{t)). 



(17) 



(18) 



This result holds in complete generality for generic Hamiltonian and transition rates, with or without conservation of 
the order parameter. In the single spin flip dynamics 



Bi(t) = 2ai(t)w([a] 



(19) 



where w([a] — > [cj]) is the transition rate between two configurations differing for the value of the spin on the site 
i. In the simulations time is discretized by single updates which, measuring time in montecarlo steps, occur on the 
microscopic time scale e = 1/N. The best numerical approximation to the limit As — > in Eq. (|17J) is obtained by 
taking As = e, which gives 



TR{t, s) = 



1 



C(t,s 



C(t,s) 



(<7i(f - e)Bi(s)) 



(20) 



In addition to R(t, s), in the following we will also be interested in the general IRF defined by 



/i(t, [t, t w }) 



R(t, s)ds 



(21) 



which corresponds to the application of the perturbation between the times t w and t < t. For this quantity, Eq. 120|) 
is replaced by 



1 e * 

Tn(t, [t,t w ]) = - [C(t,i) - C(t,t w )] - € - M - e)Bi(s)), 



(22) 



where the sum X) s =t ^ s over the discrete times in the interval [t w ,t]. 

Despite the advantages of the method, the computation of the impulsive response R(t, s) remains a very demanding 
numerical task. In order to improve the signal to noise ratio, therefore, we have found convenient to consider, instead 
of R(t,s), the quantity fi(t, [s + 6, s]) with 5/s <C 1. Assuming for R(t,s) the scaling form (J5J, to first order in S/s 
one has 



H{t, [s + S,s]) =R(t,s) 



S 2 

-h(x) + .. 
s 



(23) 



with f\(x) = (1/2) [(1 + a) + x dln /J- x ^ ]. For 5 = 1, fj,(t, [s + 5, s]) differs from R{t, s) by a correction of order 1/s, which 
we have checked to be always negligible in our simulations. Therefore, in the following, whenever results for R(t, s) 
will be presented, it is understood that the data are obtained with this procedure, if no other specification is made. 

In the following we shall use the above formulas for the d = 2 Ising model with nearest neighbors interaction and 
evolving with Glauber dynamics. 
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B. The aging contribution 

In this Section we discuss an auxiliary dynamics, referred to as no- bulk-flip (NBF), which is used 03,0111 
to isolate the aging part R ag {t, s) of the response function appearing in Eq. (J2J. 

We introduce a classification of the degrees of freedom by making the distinction between bulk and interface spins. 
A spin <Tj is regarded as belonging to the bulk of a domain if it is aligned with all its nearest neighbors. In the 
NBF algorithm bulk spins cannot flip. With this rule the interior of domains orders rapidly and all what is left 
is interface dynamics. Since in coarsening kinetics the aging contribution comes exclusively from the boundary of 
domains, measuring quantities in a simulation with the NBF rule yields the aging behavior. 

In order to illustrate this idea, let us consider the autocorrelation function. Denoting with a^ B (t) the value of the 
spin evolving with NBF dynamics, the quantity 

C NB (t,s)=M*{ar(t)crr(s)) (24) 

where M is the equilibrium spontaneous magnetization at the temperature T, is expected to coincide with C ag (t, s). 
The M 2 factor in front of the definition (|24|l is needed recalling [4j,|5|, |6j that C ag (t,s) falls from M 2 (the Edwards- 
Anderson order parameter) to zero. According to the discussion of Sec. C st (t — s) coincides with the equilibrium 
correlation C eq (t — s), which decays from 1 — M 2 to zero. In FigQ we have made the comparison between the 
stationary parts of C(t,s) computed in two different ways. In the first case we measure directly C eq (t — s) in the 
equilibrium state prepared at the temperature T/J = 2.2 (M 2 ~ 0.616). We recall that the critical temperature 
is given by Tc/J = 2.26918, where J is the nearest neighbors coupling constant. The second prescription is of 
calculating the autocorrelation function C NB (t,s) in the quench from an initial disordered state, corresponding to 
an infinite temperature, to the same final temperature T with the NBF rule, and then subtracting it from C(t, s), 
computed with the full dynamics. If the NBF algorithm correctly isolates the aging part of the dynamics what is left 
is the stationary term, C(t, s) — C NB (t, s) = C eq {t — s). However, one does not expect this to be fulfilled at all times. 
The reason is that a sharp separation between two independent components, bulk and interface, applies only when 
L(t) is sufficiently large |4j,|30j. Actually, the larger the domains are, the larger is the average distance from a bulk 
spin to the nearest interface, so that they are more effectively decoupled. The data of Fig. show a convergence of 
C(t, s) — C NB (t, s) toward C eq {t — s) as s increases, in agreement with the previous discussion. This implies that, for 
large s, the NBF algorithm does to isolate the aging part. 

One arrives at the same conclusion by considering the response function. The stationary part R eq (t — s) decays to 
zero on a characteristic microscopic time. Then, by using s and x as independent variables (which we will do from now 
on) and denoting with R(x, s) the response function in terms of these variables, in the limit of large s the stationary 
part gives a contribution only for x ~ 1. Therefore, for sufficiently large s and x > 1, one has R(x, s) ~ R ag (x, s). If 
the NBF algorithm isolates the aging contribution, the response R NB (x, s) obtained with the NBF rule and R(x, s), 
measured using the usual dynamics, should coincide for large s. In order to check this, we have computed R(x, s) 
and R NB (x, s) for a system of 1000 2 spins, quenched from the initial disordered state, corresponding to an infinite 
temperature, to the final temperature T/J — 1.5 (M 2 ~ 0.9732). We have considered several values of s ranging from 
s = 100 to s = 1600 and times t > s up to 5000. More precisely, we have used (and we shall use in the following) 
values of s generated from s n = 100+Int(1.5 n ) with n ranging from 1 to 18. The range of times considered belongs to 
the scaling regime of the system. The results in Fig. [5] show that R(x, s) and R NB (x, s) coincide within the statistical 
errors. Small differences can be detected only for the smallest values of s. This is expected since, as anticipated, the 
agreement improves with increasing s. 

In conclusion, the analysis of both the autocorrelation and the response function demonstrates the reliability of the 
NBF algorithm. Moreover, this algorithm has the advantage of speeding up considerably the simulations, since only 
the fraction of spins on the interfaces must be updated. All the results presented in the following have been obtained 
with this technique. 

IV. SCALING OF R ag (x,s) 

The first task is to check the scaling properties of R ag (x, s). If Eq. JSJ was obeyed the curves for s 1+a R ag (x, s), 
obtained for different values of s and a suitable a, should collapse on a single master curve. However, a rough 
inspection of Fig. [21 already shows that this is not the case, since the various curves cannot be superimposed by a 
simple vertical translation. We make this more precise in Fig. |3 where one can see clearly that there is no collapse 
for neither one of the two values of a proposed in Eq. I|15[l . For a = 1/4 the collapse is rather good for small x, but 
gets worst with increasing x. For a = 1/2 the collapse is bad everywhere, except for x ~ 4. 
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A. The effective exponent (Af(x, s) and the value of a 

In order to make these considerations quantitative we introduce the effective exponent defined by 

(25) 



1 + a R (x s) - - ^gg£^lf) 

Numerically, for a chosen value of x, a R ^(x, s) is given by the local slope of the plot of \nR ag (x, s) against Ins in a 
selected interval I s around s. This interval must be chosen small with respect to the range of s over which a R ^{x 1 s) 
varies appreciably. However, the smaller the interval I s the more noisy gets a^-j(x, s). On the basis of the data 
available from the simulations, the best compromise between the needs of having a local quantity and of lowering the 
noise has been reached by taking I s spanning over four consecutive values of s • The result is displayed in Fig. 
where symbols with error bars represent the numerical values of a R ^{x, s), obtained for the three different I s with the 
four values of s indicated in the legend. The continuous curves have been obtained from a fitting procedure, which 
will be discussed in Sec. IIV (Jl If Eq. (JSJ did hold, we should have found a flat effective exponent, i.e. a R ff (x, s) 



a 



independent of x and s. Instead, the data for a R ^{x, s) show the following features: 

1. Fixed s 

For fixed s, there is a strong dependence on x, revealing that Eq. (JSJ) is not obeyed. Furthermore, the curves 
display a discontinous behavior at x = 1. The equal times value at x = 1, ranging from a^-y(l, s) = —0.55 to 
a^j(l, s) — —0.53, depending on the I s considered, is separated by a jump from the smoothly increasing curve 
starting around afjy(l + , s) = 0.30. By x = 1 + we denote the smallest value of x > 1 used in the simulations. 
For larger values of x, a R jj:(x, s) keeps growing continuously. The longest set of x data, corresponding to I s with 
the smallest values of s, shows the possible saturation to an asymptotic finite value af^(oc,s) > a,f^(l + ,s). 
This behavior of a R fj explains immediately why the attempts to collapse the data with a fixed value of a, as in 
Fig. 01 do fail except in a restricted range of x. 

2. Fixed x 

The size of the error bars makes it difficult to detect an s dependence for fixed x, except in the region of short 
time difference x < 2, where with the first two I s the decrease of a R ^{x, s) upon increasing s exceeds the error 
(Fig. El- For larger values of x, error bars overlap and no statement can be made. 

Therefore, Eq. © is not obeyed and we make the assumption that the deviation is due to a correction to scaling of 
the form 

R ag (x, s) = R x {x, s) + R 2 (x, s) = s-^f{x) + s-^h(x) (26) 
where necessarily c > a, if the new term has to be subdominant. Having made this assumption, from Eq. (|25(l follows 



l + a% f (x,a) = (l + o) 
where 



1 + S- (c - a )K(x) 



(27) 



k{x) = h(x)/f(x). (28) 

Let us now pause to explore the consequences of this formula. From Eq. 1)27(1 follows that a^-j(x, s) — > a from above 
as s — > oo, for any fixed x. As pointed out previously, for those values of x where the size of errors is small enough, 
the decrease of a R j^{x, s) with increasing s is confirmed by the simulations (Fig.[SJ). Therefore, the smallest measured 
value of a^jj:(x, s) overestimates a. Excluding the value at x — 1, for the reasons which will be explained in Sec. IIV 61 
we have 

a < af//(l + ) =0.32 ±0.01 (29) 

where for a^-j(l + ) we have taken the value at x = 1 + , corresponding to the intermediate I s set. Comparing with 
Eq. I|15[) . we find that this result is compatible with Eq. I|13l) and rules out Eq. |Q . This conclusion has been 
reached with the sole hypothesis that a correction to scaling term needs to be taken into account. We have made no 
assumptions neither on the form of h(x) nor on the value of c, except for the obvious requirement c > a. In a short 
while we shall refine considerably the above estimate of a, obtaining a value much more close to the 1/4 predicted by 
Eq. O- 
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B. Analysis of Ri(x,s) and R,2(x,s) 

Let us go further with the analysis of Eq. 127|) . identifying the properties of the scaling functions f(x) and h(x), 
which must be obeyed in order to reproduce the observed features of f (x, s). Notice that R ag (x,s), being a 
response function, must vanish for large x. This requires that also f(x) and h(x) must vanish for large x. Therefore, 
the saturation to the finite asymptotic value a^y(oo, s) > a^y(l + , s), which is suggested by Fig. 01 can occur only if 
k(x) diverges for large x, that is, if f(x) decays faster than h(x). This has interesting consequences. First of all we 
get 

af ff (oo,s) = c (30) 

independent of s. Then, the correction to scaling contribution R2(x, s) is subdominant for fixed x and large s, but 
becomes dominant for fixed s and large x. Which of the two contributions Ri(x,s) and i?2(x, s) is dominant and 
which is subdominant, depends on the choice of which variable is kept fixed and which is let to grow. More precisely, 
the condition 

R%{x, s) k(x) 

defines the crossover curve x(s) which divides the (s,x) plane (FigEJ in the two regions Ti and T2, below and above 
x(s), where either Ri(x, s) or R2(x, s) is dominant. Therefore, if in the simulations one could reach values of s and t 
so large to have x(s) ^S> 1, together with a range of x extending well beyond x(s), one should observe a neat crossover 
from a^jj:(x, s) — a, in a wide interval 1 < x <C x(s) within T\, to the large- a; behavior a^y(oo, s) = c, after entering 

T2- This is visualized in the inset of Fig.^J displaying the behavior of a^^{x, s) obtained analytically from Eq. I|27|l . 
with the forms of f(x) and h(x) which will be introduced in Eqs. I|32[l and (|45|l . 
The next step is to focus, separately, on the two regions Ti and T2- 



1. Y\ region 



We make the assumption that, with the range of s reached in the simulations and x ~ 1, we are exploring the 
lower boundary of Ti, just above the s axis, where R2(x,s) is negligible with respect to Ri(x, s) (see FigEJ. The 
consistency of this hypothesis will be checked a posteriori. Let us, then, concentrate on the structure of R±(x, s). On 
the basis of existing analytical and numerical results, as discussed in detail in Ref. 01 > the scaling function in Eq. (|5j) 
is expected to be of the general form 



f(x,t /s) = A 



(x-l + t /s) a 



(32) 



where there appears the dependence on a microscopic time to, regularizing f(x, to/s) at x = 1. This extra dependence 
on to/s, for s large enough is negligible if x > 1, but becomes crucial at x = 1, where it is at the root of the observed 
discontinuity in the effective exponent. In order to see this, let us replace Eq. I|27|) with the more precise form 



■af ff (x,s) 



(l+o) 



l + 1 ^f( x ,v) + K (x,v) 

1 + S'^-^Kix^) 



(33) 



where v = to/s, n(x,v) = h(x)/ f(x, v) and the new term 



f(x, v) = v 



comes from to/s in Eq. (|32|) . A simple computation yields 



dyf(x,v) 

f(x,v) 



(34) 



_ ( -a(l - 2^1) for i-Kii 



— a— for i-l>!) 

x—l 



(35) 



which shows that this term modifies the effective exponent only at x ~ 1, while k(x,v) keeps the same properties of 
k(x), being finite for any finite x, including x = 1, and diverging for x — > oo. As a result, for s large but finite, we 
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have 

{a — a for x = 1 
a for 1 < x < x{s) (36) 
c for x 3> a?(s). 

In the simulations, the behavior at x = 1 can be studied with great precision. Considering the definition R ag {\, s) = 
limAs^o Rag{^ + As/s, s), we have computed i? ag (l, s) as i? ag (l +e/s,s), using Eq. (|20|l . This quantity is very easy to 
simulate, since good statistics can be obtained with a modest numerical effort. The results are shown in Fig.0 where, 
for completeness, we have included also the case d = 3. The data show a neat algebraic decay, consistently with the 
assumption that i?2(l, s) is negligible (if the correction was present, we should have found a crossover between two 
different power laws). Relating the slope to the top line of Eq. (|3^|) . we find 

_ f 0.473 ± 0.001 for d = 2 , . 

+ a a ~ \ 0.477 ± 0.001 for d = 3. ^ > 

These numbers need a comment. For such a short time difference x ~ -R an d the Barrat conjecture leading to Eq. i|12|) 
is correct, i.e. on the very short time regime the response function just mirrors the density of defects. Indeed, we 
have measured independently the density of defects p(s) ~ L~ 1 (s) ~ s -1 / 2 , finding 

, _ / 0.474 ± 0.001 for d = 2 , . 

/Z _ \ 0.476 ±0.001 for d = 3. ^ } 

Therefore, the comparison of Eqs. (|37|l and (|38|) leads to the identification 

l/z = l + a-a (39) 

which is, in fact, what one finds in the exact solution of the d — 1 Ising model [T2I IT3I ll7T| and in the approximate 
analytical results at higher dimensionality [l4l ITU fill . Analytical theories based on local scale invariance [23, instead, 
yield a = a + 1 , in disagreement with our data and with the aforementioned analytical results |32j . 

Once the negligibility of R,2{x, s) in the x ~ 1 region is established, the road to the direct measurement of a is open. 
This is done rewriting Eqs. JSJ and (|3"2"|l . for (t — s)/s <C 1, in the form 

s 1+a - a R 1 (x,s) = A{t- s + t )- a (40) 

which predicts, using Eq. 139II . that in the plot of s 1 ^ z R(x, s) versus t — s, the curves generated for different s ought 
to collapse, as long as (t — s)/s <C 1 is satisfied. We have made such a plot (FiglBJl of the data from the simulations 
and using the value of z in Eq. I|38|l we have found very good collapse for (t — s) < 70. From the slope of the curve 
in the region where there is collapse, we obtain a — 0.800 ± 0.005, from which, via Eqs. I|39|l and (|38|l . we get 

a = 0.273 ±0.006 (41) 

which is close enough to 1/4 to lend strong support to Eq. 1|13|) . We emphasise that this result, which has been 
obtained from the direct measurement of R ag (t,s), reproduces the previous numerical results obtained for a x from 
zfc El El 111 11. This is a nice check on the validity of Eq. (|14fl , which disproves the recent claim |26j that the 
exponent a x is unrelated to a. From the same set of data we have obtained also A = 0.153 ± 0.002J~ 1 . 



2. T2 region 

Clearly, to enter with the simulations so deep into T2 that Ri(x,s) can be neglected, is much more difficult. 
Nonetheless, the following considerations are in order. Eq. JHl requires that, for a given s, R ag (x,s) decays as 
^-•Vr/z in this sector. Therefore, the negligibility of Ri(x, s) leads to the remarkable consequence that this power law 
behavior is a property of the correction R%(x,s) and not of the leading term Ri(x, s). The implication is that the 
large x behaviors of h(x) and f(x) are given by 

h{x) - x- Xr/z (42) 
f(x) - x-^ a+ ^ (43) 

with 

a + P> \ R /z. (44) 

This cannot be accounted for by local scale invariance [2^|, which assumes that the scaling function f(x) of the 
dominant term decays like x~ Xr / z for x 3> 1. 
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C. Fitting R ag (x, s) 

Let us recapitulate what we have done so far, specifying the assumptions that we have introduced: 

• the failure of simple scaling, displayed in Fig. |3 requires to introduce a correction term. We have assumed this 
to have the form of the correction to scaling appearing in Eq. (|26|l . with c > a. This is enough to derive the 
inequality of Eq. ((29(1 . which solves the problem of the exponent a formulated in Sec. ITD 

• From the two requirements, a^(l + ,s) < eAf(oo,s) < oo and Um a; _K X , R ag (x, s) = 0, follows a^y(oo,s) = c 
and linxj^oo k(x) = oo. This, together with the assumption that f(x) is of the form ((32(1 . yields the value 141|) 
of a and the set of Eqs. ((42143144(1 regulating the large x behavior of the scaling functions. 

In order to complete the functional form of R ag {x, s), we need an additional assumption on R.2(x, s), since the power 
law behavior 1(42(1 covers only the large x behavior. However, for R.2(x, s) there are no analytical results to rely on. 
Short of any other hint, we make the simplest possible ansatz 

h(x) = Bx~ Xr/z (45) 

which continues the power law behavior 1(42(1 also into the short time region x ~ 1. Then, we arrive to the explicit 
analytical expression 

R ag (x, s) = As-^ + Bs-( 1+C h- X «/ Z (46) 

[x — 1 + v) a 

for R ag {x,s) in the al — 2 Ising model. We do not expect this to be quantitatively exact, rather with this formula 
we aim to capture the gross features of the qualitative behavior. In any case, the evaluation of the performance of 
this formula has to be made a posteriori. The program is i) to fit the remaining free parameters in Eq. ((46(1 from the 
data, ii) to check how the predicted values for R ag (x, s) compare with the numerical ones and iii) to extrapolate to 
the region of x and s that have not been reached with the simulations (see Fig[SJ). 

Our data do not allow a direct precise determination of Xr because, as Fig. |2 shows, the slope of R ag (x,s) bears 
a weak dependence on x even for the largest values of x reached in the simulations. Therefore, we take Xr = A and 
A = 1.25 [23] which is consistent, within errors, with the data. For the other parameters to,(3,B,c we have not found 
a reliable direct method to measure them. We have used a four-parameter fitting routine obtaining B — 0.47J -1 , 
c = 0.81, to = 0.01, (3 = 1.1. With these parameters, and (a, a, A) determined in the previous section llV Bl Eq. 1(46|1 
can be plotted. The comparison with the numerical data (Fig.0) is quite good in the whole range of s and x. This 
provides the a posteriori support for the validity of the procedure. 

Once f(x,v) and h(x) are given, from Eq. 1(33(1 we get the analytical form of a^^(a;,s), which can be compared 
with the data of Figs. 0] and [S] However, for a meaningful comparison, rather than making a straightforward plot of 
Eq. (|33|l . we have extracted the effective exponent from data generated from Eq. 1(4611 . following the same procedure 
described in Sec. lIV K\ Namely, we have computed R ag (x, s), from Eq. 1)46(1 . for the same values of s and t considered 
in the simulations and we have computed a^^{x : s) as the local slope over the same intervals I s of four values of s 
appearing in Fig. 0] The result is displayed in the same Figs. 0] and [5J where, although a discrete set, these values 
of a^^(cc, s) have been connected by continous lines, in order to ease the reading of the Figures. The comparison 
shows that the behavior of a^^{x, s) from Eq. 1)46(1 reproduces the numerical data within errors. However, the size of 
the error bars prevent from making claims other than qualitative. Namely, the agreement displayed in Figs.^andUJ 
between the different ways of computing a^^{x, s), makes us confident that with Eq. ((46(1 we have captured the basic 
mechanism underlying the behavior of R ag (x,s) in the d — 2 Ising model. For a detailed quantitative comparison, 
data much more precise than those presented here are needed. 

With the aim of completing the qualitative scenario, we have used Eq. l(33|) to plot a^ff(x, s), in the inset of Fig. 0] 
over ranges of s and x that cannot be accessed in the simulations. The global behavior of af^(a;, s), obtained in this 
way, displays the typical shape of an effective exponent with the crossover taking place about x(s). For s sufficiently 
large, a^^x, s) remains equal to a in the range 1 < x <C x(s), which can be enlarged at will by pushing s to larger 
values. The crossover around x(s) shows the rise of a^^{x, s) toward c. It is a very smooth crossover, since it takes 
roughly three decades for a^j(x, s) to switch from a to c. The growth of x(s) with s is rather slow and a value as 
large as s ~ 10 4 is needed in order to have a reasonably flat behavior a^-^(x, s) ~ a for x < x(s). This explains why 
the observation of a^^(x, s) ~ a over a sizable x interval is out of reach in the simulations. 
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V. INTEGRATED RESPONSE FUNCTIONS 



In Eq. (|21|l we have introduced the general form of the IRF, which contains ZFC J7J) as a particular case. The aim 
of this section is to investigate the impact on the IRF of the structure (|26[) of R ag (x, s). In particular, we shall have 
to understand why the existence of such a strong correction to scaling in R ag (x, s) does not surface at the level of 
ZFC, where the simple scaling form © accounts very well for the data We shall not discuss the behavior of 

TRM because, as mentioned above and explained in detail in |0, , this quantity is the most unfavorable choice for 
the analysis of the scaling properties. 

Inserting Eq. I)26|l into the definition H21|) and keeping track of the dependence on to, we get 

v ag (y,t w ) = txi{y,t w ) + v 2 (y,tw) = t~ a F( y ,v) + t- c H( y ) (47) 

where 

f u 

F(y,v) = y- a z-( 1+a U(z,v/y)dz (48) 



with y = t/t w , u = t/t, v = to/t w and 



H{y) = y- C f 
Ji/y 



(z,v/y) = A 



,-(!+<=) 



h{l/ z)dz 



v a+f3 



(1 - z + v/y) a 



(49) 



(50) 



comes from Eq. <|32[) substituting x with l/z. For simplicity, we have omitted to indicate explicitely the dependence 
on the upper limit of integration u. The effective exponent 



ae ff {y,t w ) = - 



din /j, ag (y,t v 



din t u 



-(c-o) 



K(y,v) 



i+t; 



-(c-o) 



K{y,v) 



has the same structure of Eq. , with 



and 



F{y,v) 



K(y,v) 



d v F(y,v) 
' F(y,v) 

H(y) 

' F(y,v) 



(51) 



(52) 



(53) 



being the analogues of f(x,v) and k(x, v). 

There are as many IRF as there are ways of choosing t. Here, we restrict the attention to the two cases with 
t/t w = q and t/t = p, where q and p are fixed numbers. It will be convenient to introduce the notation aj ag and ir ag 
for these two particular IRF. The distinction between the two cases enters the above formulas in the upper limit of 
integration in Eqs. I|48() and (|49|l . where 



q/y 



for 



for 



'ag- 



(54) 



We shall now establish the gross features of a^j:j:(y,t w ), as y is varied with fixed t Wl using the properties of F(y,v) 
and K(y, v) derived in the Appendix. The behavior of these functions depends on the sign of a — 1. With the value 
of a = 0.80 obtained from the fit of Eq. (|40|l . we shall make use only of the results for a < 1, leaving the general 
discussion to the Appendix. 
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A. LU ag 

With y > t/t w — q, where q > 1 is some fixed number, y takes values in [q, oo). From Eqs. and l(7T|) in the 
Appendix follows 

while, from Eqs. (|72[) and (|76|) follows that K UJ (y,v) is finite for finite y and diverges for y — > oo. Hence, the affective 
exponent 

behaves like af/j in Eq. I|36|) . apart from the absence of any discontinuity at the minimum value of y, since Eq. I|55|l 
shows that F u (y,v) vanishes for all y as t w — > oo. This is displayed in FigGO illustrating the behavior of a%fAy,t w ) 
obtained from the simulations with t w in the range [200,400]. 

B. TY ag 

If we keep fixed the ratio t/t = p, with p < 1, y takes values in [1/p, oo). In this case we obtain an IRF in which the 
correction to scaling does not produce the crossover in the effective exponent from a to c, as observed in the previous 
case and in R ag . Rather, a^jfiy, t w ) starts from a and ends up again to a, as y grows, with a possible discontinuity at 
y = 1/p, as will be explained below. The difference comes from Eq. I|82ll in the Appendix, showing that now K 1T (y^v) 
vanishes as y — * oo. After taking into account that also F n (y, v) vanishes as y — > oo (see Eq. (|fi3*|) in the Appendix), 
this implies 

lim al ff (y,t w ) = a. (57) 

y — >oo J J 

Around the minimum value y = 1 /p, the behavior is different for p < 1 and p = 1 . In the former case no singularity 
develops, since from Eq. (|8t)(l follows that F v {y = l/p,v) ~ v. Instead in the latter case, which corresponds to ZFC 
(i.e. ir ag — Xag for p = 1), from Eq. Ij84[l we have that at y = 1 there is a discontinuity in the effective exponent with 

a* ff (l,v,t w ) = a - a (58) 

exactly as in Eq. ijjffijl at x = 1. 

The overall behavior of a^j(y, t w ) is displayed in FielTniobtained by plotting Eq. l|5"T|l for fixed t w . The two curves, 
corresponding to t w = 200 and t w — 2000, display a very fast rise from a — a — —0.53 at y — 1 followed by a very 
slow approach to the asymptotic value l|57() . The absence of the crossover is the most prominent qualitative difference 
with respect to the inset of Fig 01 which explains why the correction to scaling plays a minor role in the analysis of 
ZFC data. This very mild increase of a^^(y,t w ) in the range 1 < y < 10 is observed also in the simulations (notice 
the vertical scale in Fig lll|) . 

Furthermore, notice that in Fig llOl the asymptotic value is reached from above for t w = 200 and from below for 
t w = 2000. The two curves cross each other, due to the interplay of the relative weights of F x (y,v) and K x (y,v) as 
t w is varied, the former being a negative quantity and the latter a positive one. Now, the important point is that this 
particular feature can be resolved also in the data from the simulations (Fig. Ill|) . showing that Eq. (|46|l does, indeed, 
account for the observed behavior of ZFC. 

VI. CONCLUSIONS 

In this paper we have studied numerically the linear response function of the d = 2 Ising model quenched below 
the critical temperature. The data for R ag show that the simple scaling form {5J, usually believed to hold in the 
aging regime, is not obeyed. We have attributed the deviation from simple scaling to the existence of a correction 
to scaling R.2{x, s), which, in order to explain the observed behavior of the effective exponent, although subdominant 
for fixed x and growing s, must become dominant for fixed s and large x. Then, focusing on the time regime where 
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R,2(x,s) is negligible, we have been able to measure a with good precision obtaining a = 0.273 ± 0.006. This solves 
the problem of the exponent a formulated in Sec. [n] and is in agreement with previous numerical results from the 
ZFC. Furthermore, this value is well consistent with Eq. (|13[) and confirms the dependence of a on dimensionality, as 
found in all cases where analytical results are available. 

Awareness of the existence of this correction to scaling is of fundamental importance when analyzing numerical 
data, since an interpretation which does not take into account the strong dependence of a^j(x,s) on x may lead 
to wrong conclusions. In particular, if insisting in collapsing the data according to Eq. JSJ, one would find that the 
best data collapse obtains with an exponent whose value is somewhere between a and c, depending on the range 
of x considered in the simulations. However, the collapse cannot be satisfactory in the whole range of x, as Fig. [3| 
shows. The findings of Chatelain [l8j. perhaps, can be interpreted along these lines. Indeed, in Ref. jl8j the impulsive 
response function was computed and the best data collapse was obtained with a — 1/2 using data for rather small 
values of s. The author did explicitly notice that the collapse was not perfect and hypothesized himself the existence 
of strong corrections to scaling. 

We have also investigated the functional form of the scaling function entering i?i(x, s), and we have found that 
Eq. (j32J) with a = a + 1/z, which was proposed in [l7l | as a generalization of the known analytical results, compares 
quite well with the numerical data. Next, by making the simplest possible ansatz for i?2(x, s), we have obtained an 
analytical expression for the full R ag ■ This is quite useful for the investigation of the overall qualitative behavior and 
compares satisfactorily well with the data in the time regions reached by the simulations. 

Finally, we have considered the retrival of the properties of the impulsive response function from the IRF. As 
discussed in 01 > this may be an issue to treat with care when dealing with ZFC, due to the presence of to, which acts 
as a dangerous irrelevant variable for a > 1 or, equivalently, for d > djj. However, in this paper we are not concerned 
with this aspect of the problem, since d = 2 is below djj = 3 and we have a = 0.80 < 1. Nonetheless, recovering the 
scaling exponent may still be complicated, due to the presence of the correction i?2(£, s). In fact, we have found that 
the relative importance of this term strongly depends on the kind of response function considered. In particular, it is 
maximal in R ag (t,s), while it is almost negligible for Xag(t,s). This implies that i) by weighting differently Ri(t,s) 
and i?2(i, s), different response functions may behave very differently and ii) the best suited function in order to 
weaken the correction term and to access the asymptotic properties is Xag(t, s). 

The rich pattern of behaviors uncovered in this paper shows that there remains much to be understood of the 
response function in slow relaxation phenomena, even in the relatively simple case of coarsening systems. Presently, 
we do not know how general is the behavior that we have found in d = 2. It would be interesting to perform simulations 
in higher dimension to address this point, at least numerically. 



VII. APPENDIX 



In order to analyse the behavior of 

a^ ff (y,t w ) = a 



l + lF(y,v) + H~ {c ' a >K(y,v) 
l + t- {c - a) K(y,v) 



(59) 



it is first necessary to study the properties of F(y, v) and K(y, v). Rewriting Eq. I|48|) as 

F(y,v)=Ay- a I(a) (60) 

with 

ru p +a -(i +a ) 

1(a) = / dz- — -p- (61) 

Ji/y {l~z + v/y) a 

we have 

d v F{y, v) = -Aay- [1+a) I(a + 1) (62) 

and 

F(y,v) = -a[-) \ 63 

Next, carrying out the integration in l|49|l . we can write 

K(y, v) = B^- {u x «/ z - c - y-^l*] (64) 
1(a) 
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where B = , A B , r . 

The integral 1(a) is a finite function of v/y everywhere, except at v/y = 0, where a singularity develops if a > 1 
and u — ► 1. In order to see this, let us separate in /(a) the contribution coming from the neighborhood of the upper 
limit of integration, by writing 1(a) = Iq + I u , which corresponds to the partition of the domain of integration 

U — € ru 

+ / (65) 
where e is a small number. The first contribution Iq is always finite, while 

Iu = L dz d- z+ v/ y r (66) 



/3+a-(l+a) 



-[(1 - u + e + tVy) 1 -" — (1 — t* + w/y) 



1-Q1 



1 — a 

is finite for u < 1, but diverges for u = 1 and t>/y — * oo, if a > 1, yielding 

!< oo for a < 1 
"M^T for a=l (67) 
^(f/y) 1 -" for a>l. 

Going through similar steps, one can see that also I (a + 1) is always finite for u < 1, while for u = 1 and small i>/y 

7(a + l)~i( W /y)- Q . (68) 
a 

Therefore, from Eq. I|63(l and for u < 1 we have 

F(y,w)~-(«/y) (69) 

while for u = 1 

-(v/y) 1_Q for a < 1 

«)H R^) for a = 1 ( 7 °) 

1 — a for a > 1. 

The above results hold in general. In the following we continue separately the discussion of ui ag and 7r Qg . 

A. LO ag 

With 1 < q < y and u = q/y, the integration domain [l/y,q/y] in Eqs. (|48(l and l|49(l includes u = 1 only when y 
takes the smallest possible value y = q. Therefore, from Eq. H7U|) we obtain 

-(v/?) 1 -" for a < 1 

*Ug,«)H rW for a = 1 ( 71 ) 

1 — a for a > 1 



while from Eqs. JMIl and (|BTjl 



< oo for a < 1 

4 -R^y for a = 1 (72) 

(w/g)"" 1 for a > 1. 



Inserting these results into Eq. 159(1 . for the effective exponent at y = q we get 

a for a < 1 
a + 1 — a = 1/z for a > 1 



H ™ *») = i a + 1 = 1°/^ 1 for r> > 1 . ( 73 ) 
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Instead, for y > q, since F(y,v) obeys Eq. i(6"5|) and K LJ (y,v) is finite, we get 



lim a% ff (y,t a ) = a (74) 

t w — *oo J J 

independently from the sign of a — 1. However, care must be used when y — > oo, because the domain of integration 
in Eq. (|61|l shrinks to zero and 1(a) vanishes like 

1(a) ~ y a ~^+ a ) (75) 

yielding 

iUy, v)~/ +a - Wz - (76) 

Therefore, from Eq. (|44|l follows that K UJ (y,v) diverges for large y and, in turn, this gives the non commutativity of 
the limits lim^oo and limt M _ >oc ; since keeping t w fixed 

lim a" * Ay, tw) = c (77) 

y — >oo 

while with y > q fixed 

lim a" ff (y,t w ) = a. (78) 



B. n ag 

The integration domain [l/y,p] in Eqs. I|48|l and l|49|) is given by [l/y,p], with y > I /p. Therefore, the upper limit 
of integration u = 1 is not reached if p < 1, while it is reached if p = 1. Furthermore, as y hits the lowest value 
y = 1 /p, the integration domain shrinks to zero and to first order in (y — 1 jp) 

1(a) = p 2+p+a - {1+a \l ~ p(l - v)]- a (y ~ l/p) (79) 

from which follows 

lim F«(y, v) = -a PV . (80) 

y^l/p 1-p+pv 



Notice, also, that from Eq. (|rjl|) we get 



KAy, «) = t£t [/ R/Z - C y a - C - y- W«] (81) 
/(a) 

where both the exponents of y in the right hand side are negative. Hence, for fixed t w 

YimK„(y,v) = Q (82) 

y^oo 

contrary to what happens in the case of w ag , where the limit for y — > oo produces a divergence. Let us proceed by 
considering separately p < 1 and p = 1. 

p < 1 

In this case u < 1 for all y and /(a) does not develop any singularity. From Eqs. I|69(l . 1)80(1 and (|81f) then follows 

lim a^ ff (y,t w ) = a (83) 

for all y. 

p = l 



From Eq. lf%U|) . at y = 1 we have 



—a 



(84) 
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while for y > 1 Eq. (|70|l applies. Switching to the behavior of K^(y,v), notice that inserting Eq. I|67(l into Eq. I|81|) 
we obtain 

y a-c _ y a-X R /z for a <1 

K„(y,v)^( - va ~°ar^ for a = l (85) 



ln(v/q) 

( y a-c _ ya-Xx/z^/^a-l for a> 1. 



Using these results, for fixed y we get 



while for fixed t, : 



a — a — 1/z — 1 for y = 1 
lira a^j:(y,t w ) = ( a f° r 2/ > 1 an d a < 1 (86) 
x a + 1 — a = 1/z for y > 1 and a > 1 



,'i n i "'<'/'<"'• '"'> i o+i? r a = l% <1 for a>l. (87) 



y^oc 



Finally, let us recall that from Eqs. i|13|) and l|39fl we have 



1 1 d-du faa _ 

a-l = -- — (88) 

zdu — at 



which shows how the sign of a — 1 changes with the dimensionality. 
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t-s 

FIG. 1: Plot of the difference C(t, s) - C NB (t, s) for s = 6 ■ I0 2 , 2 • fO 3 , 10 4 , 5 • fO 4 (from bottom to top) showing the convergence 
toward C eq (t — s) ( bold curve) as s increases. 




FIG. 2: R ag (x,s) (continuous lines) obtained with the NBF algorithm and R(x,s) (squares) obtained with the full quench 
dynamics are plotted for s = 111, 229, 537, 1085, 1577 from top to bottom. Circles represent the fit according to Eq. 14611 . 




FIG. 3: Failure of collapse of the curves of Fig.|2]witli a = 1/4 (lower set) and a = 1/2 (upper set). 
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FIG. 4: Symbols with error bars represent afjf(x,s) obtained from the simulation data, as local slopes of In R(x,s) against 
Ins over the four values of s specified in the legend. Continuous curves are the plot of c&Ax,8) obtained from the fitting 
formula (1461 for the same s used in the simulations. The frame is magnified in Fig. [5] The inset shows a^^{x,s) from 
formula I46H over a range s = 3 • 10 2 , 10 3 , 10 4 , 10 5 , 10 6 (from top to bottom) that cannot be reached in the simulation. 




FIG. 5: Magnification of the framed portion of Fig.|lJ showing that the decrease of aft Ax, s) with increasing s (limited to the 
first two sets I s ) exceeds the error bars for x < 2. Symbols are the same as in Fig. [I] 




FIG. 7: Plot of R ag (l, s) in d = 2 and d = 3. Straight lines are power law best fits £-°- 473 in d = 2 and i" - 477 in d = 3. 




FIG. 10: <x*ff(y,tw) plotted from Eq. 15 11 with t w = 200 (squares) and t w = 2000 (circles). The horizontal line corresponds to 
a = 0.273. E 



